function avgOut = simulator_threshold_finite(theta2,output,thresholdLevel,numSim)
% Here we simulate repeated small samples and report the average to correct
% for finite-sample biases. 
minObsPerRegime = 15;
T               = size(output.factors,2);
xHat            = [output.macros';output.factors];

% We simulate for all maturities
tmp      = output.model;
setup.nn = tmp.nn;
setup.nm = tmp.nm;
nm       = tmp.nm;
setupAll           = setup;
setupAll.matSelect = 1:max(tmp.matSelect);
if output.model.modelNum == 1
    [modelAll,errorMes]= solveQTSM(tmp.alpha,tmp.beta,tmp.psi,tmp.mu,tmp.phi,tmp.sigmaxtil,setupAll);
elseif output.model.modelNum == 2
    [modelAll,errorMes] = solveQTSM_constRsmooth(tmp.alpha,tmp.rhor,tmp.beta,tmp.psi,tmp.mu,tmp.phi,tmp.sigmaxtil,setupAll);
end
modelAll.modelNum = output.model.modelNum;
modelAll.factorSignResOn = 0;
modelAll.rhor     = output.model.rhor;

% Unfold the factor dynamics
nx = size(modelAll.gx,2);
[h0_1,h0_2,hx_1,hx_2,sigma,gama0,gamaz,gamax,sigmaz]   = unfoldTheta2_threshold_macro(theta2,nx);

%% Simulation
rng(1,'twister')

dimy           = size(modelAll.g0,1);
modelAll.macros= zeros(T,nm);
modelAll.factorSignResOn = 0;
out            = cell(numSim,1);
s = 0;
while s < numSim
    if mod(s,50) == 0
        disp(['Number of draws: s = ',num2str(s)])
    end
    xSim           = zeros(nx,T);    
    zSim           = zeros(T,1);
    ySim           = zeros(dimy,T);
    rLag           = 0;
    startIdx       = randi(T,1,1);
    xSim(:,1)      = xHat(:,startIdx);
    zSim(:,1)      = output.econAct(startIdx);
    for t=2:T
        if zSim(t-1,1) >= thresholdLevel
            xSim(:,t) = h0_1 + hx_1*xSim(:,t-1) + sigma*randn(nx,1);
        elseif zSim(t-1,1) < thresholdLevel
            xSim(:,t) = h0_2 + hx_2*xSim(:,t-1) + sigma*randn(nx,1);
        end
        zSim(t,1) = gama0 + gamaz*zSim(t-1,1) + gamax'*xSim(:,t-1)+ sigmaz*randn;
        modelAll.macros(t,:) = xSim(1:nm,t)';
        modelAll.rLag(t,1)   = rLag;
        ySim(:,t) = modelFit(modelAll,xSim(nm+1:end,t),t);
        rLag = ySim(1,t);
    end
    % Test if we have both booms and recressions in z_Boot
    regime_1       = zSim >= thresholdLevel; %a boom
    regime_2       = zSim < thresholdLevel;  %a recression
    if sum(regime_1) >= minObsPerRegime && sum(regime_2) >= minObsPerRegime
        s = s + 1;
        % The regressions
        ySimTrans = ySim';
        h_index = [1,3,6,9,12,15,18,21,24];
        endIdx  = max(output.model.matSelect);
        for i=1:length(h_index)
            shortRateOneYear  = ySimTrans(:,modelAll.matSelect==12);
            startIdx          = h_index(1,i);
            maturities        = modelAll.matSelect(1    ,startIdx:h_index(1,i):endIdx);
            yields            = ySimTrans(1:end,startIdx:h_index(1,i):endIdx);
            shortRate         = yields(:,1);
            
            % Moments for model evaluations
            tmp = ['h',num2str(h_index(1,i))];
            bootOn = 0;
            out{s}.(tmp) = momentsModelEval(yields,shortRate,shortRateOneYear,zSim,modelAll.macros(:,2),thresholdLevel,maturities,h_index(1,i),bootOn);
            %out{s}.(tmp) = rmfield(out{s}.(tmp),'dataDyieldRec');
        end
    end
end

% Computing averages 
namesPartI  = [{'meanYield'},{'meanYield_regime1'},{'meanYield_regime2'},...
                {'stdYield'},{'stdYield_regime1'},{'stdYield_regime2'}]; %{'avgdYieldsRec'}
avgOut = out{1};
for h=1:length(fieldnames(out{s}))
    for n=1:length(namesPartI)
        name = namesPartI{n};
        for s=2:numSim
            avgOut.(['h',num2str(h_index(h))]).(name) = ...
                avgOut.(['h',num2str(h_index(h))]).(name) + out{s}.(['h',num2str(h_index(h))]).(name);
        end
        avgOut.(['h',num2str(h_index(h))]).(name) = avgOut.(['h',num2str(h_index(h))]).(name)/numSim;
    end
end
% for structures
namesPartA = [{'LPYI'},{'LPYI_Regime'},{'xhrOnSlope'},{'returnsOnShortRate'},...
    {'pmiOnSlope'},{'ygapOnSlope'},{'AR1_regime'}];
for h=1:length(fieldnames(out{s}))
    for j=1:length(namesPartA)
        nameA = namesPartA{j};
        namesPartB = fieldnames(avgOut.(['h',num2str(h_index(h))]).(nameA));
        for n=1:length(namesPartB)
            nameB = namesPartB{n};
            for s=2:numSim
                avgOut.(['h',num2str(h_index(h))]).(nameA).(nameB) = ...
                    avgOut.(['h',num2str(h_index(h))]).(nameA).(nameB) ...
                    + out{s}.(['h',num2str(h_index(h))]).(nameA).(nameB);
            end
            avgOut.(['h',num2str(h_index(h))]).(nameA).(nameB) = ...
                avgOut.(['h',num2str(h_index(h))]).(nameA).(nameB)/numSim;
        end
    end
end


end